Lattice QCD study of a five-quark hadronic molecule 
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We compute the ground-state energies of a heavy-light K-A like system as a function of the 
relative distance r of the hadrons. The heavy quarks, one in each hadron, are treated as static. 
Then, the energies give rise to an adiabatic potential V a (r) which we use to study the structure 
of the five-quark system. The simulation is based on an anisotropic and asymmetric lattice with 
Wilson fermions. Energies are extracted from spectral density functions obtained with the maximum 
entropy method. Our results are meant to give qualitative insight: Using the resulting adiabatic 
potential in a Schrodinger equation produces bound state wave functions which indicate that the 
ground state of the five-quark system resembles a hadronic molecule, whereas the first excited state, 
having a very small rms radius, is probably better described as a five-quark cluster, or a pentaquark. 
We hypothesize that an all light-quark pentaquark may not exist, but in the heavy-quark sector it 
might, albeit only as an excited state. 

PACS numbers: 12.38.Gc, 12.39.Mk 
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INTRODUCTION 

Lattice QCD studies of hadron-hadron interactions are 
the gateway to nuclear physics through first principles 
1] . From a lattice simulation point of view the nucleon- 
nucleon interaction is unquestionably the most challeng- 
ing case 0], and it might not be resolved in the foresee- 
able future. However, interactions in other two-hadron 
systems are worth investigating as well, because new in- 
sights into the structural features of already discovered 
[3| or yet unknown baryon resonances may emerge. In 
particular, one may ask if some of those may be under- 
stood as hadronic molecules, like the deuteron, or if more 
compact clusters, like a pentaquark may also exist. 

We are here interested in pairs of hadrons contain- 
ing one heavy quark each. In such systems a rela- 
tive, residual, interaction is a well-defined concept. In 
the spirit of the Born-Oppenheimer approximation, the 
(slow) heavy quarks serve to define the centers of the two 
hadrons while the (fast) light quarks and gluons provide 
the physics of the interaction. Exploratory studies along 
those lines have been done in the context of meson- meson 
and baryon-baryon systems [1, 0] • 

Specifically, we here investigate a heavy-light meson- 
baryon (five-quark) hadron with the quantum numbers 
of an S-wave K-A system. The heavy-quark propagator 
is treated in the static approximation. This allows us 
to compute the total energy as a function of the relative 
distance r between the heavy quarks, viz. hadrons. A 
production scheme is illustrated in Fig. [TJ The result- 
ing adiabatic (Born-Oppenheimer) potential V a (r) then 
may be used to address the possibility of molecule-like, 
or other, structures. 

In the static approximation to the heavy-quark prop- 
agator the potential V a (r) extracted from the lattice 
equally applies to the systems K-A, D-A c , and -B-A&, 




FIG. 1: Schematic view of K-A like molecule production and 
illustration of the current scheme for extracting a hadronic 
interaction at relative distance r. Thick lines indicate heavy- 
quark propagators, and thin lines depict light quark propaga- 
tors. 



the s-flavored quarks being replaced with c and b flavors, 
respectively. Within the framework of a Schrodinger 
equation those systems can be studied using their respec- 
tive physical reduced masses, provided of course, one is 
willing to accept the limitations of the adiabatic approxi- 
mation and the non-relativistic nature of the framework. 
This is the main reason for classifying our results as qual- 
itative. 



OPERATOR CONSTRUCTION 

The lattice simulation requires two-hadron interpolat- 
ing fields. Those are constructed using as building blocks 
standard local operators for the K + and A particles 0] . 
They are placed at relative distance r and then projected 
to total momentum zero 

O a (r; t) = A- 1 / 2 £ ]T V, f - yK+mAlm ■ (1) 
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Here the normalization N is the spatial lattice volume, 
the sums are over lattice sites, and a is a Dirac spinor 
index. Then, with 

d^r;t) = Oi{r;t)^ >alJt (2) 

the time correlation function 

C(t,t ) = (O^r:t)0^s-M)) - (O t)) (6 „(S| t )) , 

(3) 

where s is the relative distance at the source, can be 
expressed in terms of fermion propagators. The flavor 
assignment K + A ~ su uds causes the separable term in 
([3]) to vanish. Writing H(xt,yto) and G(xt,yto) for the 
heavy (s) and light (u,d) quark propagators, respectively, 
one obtains 

C(*,*o) = (Ey[H(yt,y + rt)H{n + 810,7^0) 
-H(yt, rit )H{?x + st , y + ft)} x 

G(yt, r!to)[G(yt, r x t )G(y + ft, r x + st ) 
-G{yt,n + st )G(y + rt,nt )] ) . (4) 

For clarity the rather involved color and spin index struc- 
ture is not shown in ([4]). Also, translational invariance of 
the gauge field average ( ) has been used to arrive at the 
above expression, and utilizing this freedom, an arbitrary 
space site f\ was introduced to fix the source locations. 
A diagrammatic representation of ((3]) is shown in Fig. [5] 




FIG. 2: Diagrammatic representation of ((4]). Thick lines in- 
dicate heavy-quark propagators, and thin lines depict light (u 
and d) quark propagators. 

The heavy-quark propagators are employed in the 
static approximation. For (unimproved) Wilson fcrmions 
with hopping parameter k this means that the propaga- 
tor is taken in the limit k — > , resulting in 

H{xt,yt ) = 8 s ^ (2/t)*- to i(l + 74 )Wt(j?;tot), (5) 



where l4(S;tot) is the product of SU(3) link variables 
along a straight line from (xfo) to (aft) 0]. The factor 
(2K) t_t ° gives rise to only a constant global energy shift 
Am = — ln(2/«) , which we will ignore. 

The distance r — is rather special [6(] because a color 
singlet operator, as realized by |T]), can also be achieved 
by a "color twisted" version of |T]) where quarks across 
the hadrons K + and A , now at the same location, are 
combined into a color singlet. Because we do not con- 
sider color twisted operators in this work, we restrict 
ourselves to non-zero relative distance [26j]. Thus, using 
H(yt,y + ft) oc and H(f\ + sto, fitg) oc the first 
two terms in ((4]), and accordingly the top two diagrams in 
Fig. [2] vanish, and only the last two of those make a con- 
tribution to the correlation function for non-zero relative 
distance. By way of ([5]) those contributions are propor- 
tional to Sff^Sps- Introducing f*2 = ?x + s, to replace s, 
the correlation function (j4j becomes 

C(t,t Q ) = S^ f2 _ fi (H{? 1 t,r 1 t Q )H{? 2 t ,r f 2 t) (6) 

Gint, r*it )[-G(rit ) f^G^t, f 2 t Q ) 
+G(f 1 t,f 2 t )G(f 2 t,f 1 t )}) for rVO. 

As a consequence of the static approximation, the site 
sum J^y h as vanished from and thus, is unfortunately 
no longer working to improve statistics. 

The final correlator we use is extended from ([6]) to a 
K x K matrix by employing several levels k — 1 ... K 
of operator smearing. The procedure amounts to re- 
placing in {T]) all light-quark fields ip,ip with smeared 
fields ij){ k } T%j){ k } . We have used Wuppertal-style fermion 
smearing [9( along with APE-style gauge field fuzzing 
[Io| . The implementation specific to the asymmetric lat- 
tice used here is discussed below. No smearing, nor link 
variable fuzzing, was done for the heavy, static, quark 
fields in order to preserve spatial locality, i.e. the S fac- 
tor in |5j). Thus, generically, replacing O — ► 0^ k ' = 
correlator (J6j) becomes a K x K ma- 
trix 

C kk " \t,t )(r;t,t ) = (OW(f;t)6& >(F;t„)>, (7) 

with k, k' = 1 ... K and a sum over fi is understood. 
The expression for C kk (t,to) in terms of quark prop- 
agators still has the form given by ([6]), however, light- 
quark propagator elements are replaced with smeared 
ones, G — > G^ kk h Smearing and fuzzing prescriptions 
are the same at source and sink. Thus the correlation 
matrix ([7]) is hcrmitian by construction. 

We choose an asymmetric and anisotropic lattice with 
geometry L\ x L 2 x L 3 x L 4 = 8 x 8 x 32 x 16. The 
bare lattice constants, in the respective directions, satisfy 
ai = a 2 = 2a3 = 2a4. Subsequently we shall refer to 
«4 = «3 = a as the lattice constant a of the simulation, 
so the physical volume is 16a x 16a x 32a x 16a. The 
lattice has a fine resolution in 4-direction (time) , and the 
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same fine resolution in the (spatial) 3-direction where the 
lattice size is twice as long as in the 1,2-directions. The 
idea here is to place the static quark sources along the 
3-dircction and thus achieve a fine spatial resolution of 
the relative distance as well as to provide enough space 
to allow for asymptotic separation of the hadrons. 

The spatial lattice asymmetry leads us to modify the 
procedure for operator smearing. A smearing iteration of 
k steps is understood to mean k steps in the 1,2-directions 
but 2k steps in the 3-direction, so the physical smearing 
extent is the same for all spatial directions (neglecting 
renormalization effects). Mutatis mutandis, the same is 
true for gauge field fuzzing. An illustration is shown in 
Fig. [3] The implementation used here is described in 
more technical terms in reference 



Smearing direction 



Fuzzing direction i 



FIG. 3: Illustration of the prescription for quark field smear- 
ing (top) and gauge field fuzzing (bottom) on the asymmetric 
lattice. The geometry for k — 1 iteration is shown. 

The positions of the static quark sources are chosen as 
x = (5, 5, n, 3) with n = 1, 2, 3, 4, 8, 11, 13, 17. Their spa- 
tial locations along the 3-direction, see Fig.[5J allow us to 
achieve any relative distance r = la . . . 16a by choosing 
an appropriate pair of source points. Some relative dis- 
tances can be obtained more than once. In those cases ap- 
propriate averages over sets of source and sink pairs have 
been taken to compute the correlator matrix elements. 
Because periodic boundary conditions are applied in the 
spatial directions relative distances larger than r = 16a 
are redundant. 
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FIG. 4: Placement of static hadron sources along the 3- 
direction of the lattice. Periodic boundary conditions are ap- 
plied to the quark and gluon fields. All distances between 
r = la and r = 16a can be achieved, some multiple times. 

Given the approximations employed in this work, e.g. 
static heavy quarks, and a few more mentioned below, 
it would be of no advantage to employ sophisticated im- 



proved lattice actions. We have therefore used the stan- 
dard Wilson plaquette action with Wilson fermions in 
a quenched simulation Q. The gauge field couplings in 
the p-v planes and the hopping parameters in directions 
p are given by, respectively, 



aia 2 a 3 a4 



and 



(8) 



The simulation was done at (3 = 6.2 with four values 
K = 0.140, 0.136, 0.132, 0.128 of thehopping parameter 
employing a multiple mass solver 



lj. A total of 343 



gauge field configurations were used. 



MASS SPECTRUM ANALYSIS 

The time evolution of the eigenvalues of C(t, to) deter- 
mines the energy spectrum. Eigenvalues, behaving expo- 
nentially with t, may rapidly vanish into numerical noise. 
Conventional diagonalization methods do not work well 
under those circumstances. Singular value decomposition 
(SVD), on the other hand, is ideally suited to the task 
[23]. The SVD reads 



C(Mo) = U(t,t )E(t,t )vHt,t ) 



(9) 



where U(t,to) and V(t,to) are unitary in our case[27j. 
and £(i,£o) — diag(<7i(t, to) • • • <?K (^o)) contains the 
singular values satisfying ak(t,to) > 0. If C(f, <o) is 
non-degenerate and positive semi-definite then the set 
of singular values {<7fe(t, to) : k — 1 . . . K} and the set 
of eigenvalues are the same. A few more details can be 
found in 11|. For simplicity we will refer to <7fe(t,to) as 



eigenvalues. 

To extract energy levels from the eigenvalues an often 
used procedure, sometimes starting from a generalized 
eigenvalue problem [l4j], is to construct effective mass 
functions from the eigenvalues and then try to identify 
plateaus in the asymptotic time region. Typically, only 
a narrow subset of the time slices is usable, which also is 
subject to some discretion. We shall not rely on effective 
mass functions here, but rather analyze correlator eigen- 
values by means of the maximum entropy method [l5| fol- 
lowing the adaption described in [ll[ , and in [l6| . Among 
the advantages are that all available time slices may be 
used, if desired, and that excited states, if present, are 
easily revealed. 

The time dependence of the eigenvalues <Xfc(t,to) shall 
be fit with the model 



F(p\t) 



dui p(ui) 



(10) 



for some set of discrete times slices t — t\ . . . t%, where 
p(u>) is a spectral density function. The maximum en- 
tropy method (MEM) is based upon Bayesian statis- 
tics. In this context the numbers p{u>) are interpreted 
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as stochastic variables [la ]. Their most likely values are 
obtained by simulated annealing as described in [llj] . We 
refer the reader to this reference for technical details. 

At each of the 16 relative distances r we have used 
three (K = 3) operator smearing levels k = 1,2,3. 
A typical set of eigenvalues is displayed in Fig. [5] for 
distance r = la and the smallest (light) quark mass, 
K = 0.140. Only two eigenvalue correlators are shown, 
the third one is indistinguishable from zero due to nu- 
merical noise. Note that eigenvalues are separated by 
about four orders of magnitude. Thus crossing ambigui- 
The eigenvalue sets for all other 



ties are nonexistent 
distances r exhibit the same features. 



28] 




14 16 



FIG. 5: Two sets of eigenvalues ak(t,to),k = 1,2, of the 3x3 
time correlation matrix © for relative distance r — la and 
quark hopping parameter k — 0.140. The third eigenvalue, 
k — 3, is not shown, it vanishes into numerical noise. Filled 
plot symbols belong to time slices used for the MEM fits em- 
ploying the model (|10|) . the curves are the corresponding re- 
sults. Note that the source has been shifted to to = 0. 



The eigenvalue correlators were ht with the model (jlO|) 
using, consistently in all cases, time slices t = 1 ... 8 
with the source tg = being omitted, and mass range 
of < u> < 4 with discretization interval size Aw = 0.05, 
all in units of a -1 . The maximum entropy method then 
yields spectral density functions p(lo). Representative ex- 
amples are shown in Fig. [6] for r = 1,4, 8, 16, in units of 
a, and k = 0.140. The histogram lines in Fig. [6] repre- 
sent p(u>), obtained from simulating annealing (cooling) 
following the very procedure put forward in 11[ . Sixteen 
random annealing starts were used. The solid histogram 
lines represent the average of those 16 runs, the dashed 
histogram lines give the corresponding standard devia- 
tion. In most cases the spectra exhibit isolated peaks, 
say S n = {uj : u> G peak #n}. Then one may compute, 
for each peak n, the volume Z, the energy E, and the 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 




3.0 3.5 4.0 




0.0 0.5 



IS 3.0 3.5 4.0 



FIG. 6: Spectral density functions p, according to (|10[1 . for 
two-hadron operators with relative distance r = 1, 4, 8, 16, in 
units of a, thick histogram lines. Shown are eigenvalue spectra 
which contain ground state peaks. The dashed histogram 
lines indicate errors as explained in the text. The smooth 
dotted lines correspond to Gaussian fits to the discretized 
spectral densities. 
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width A, according to 

Z = J dujp(uj) (11) 

E = Z' 1 J dup{u)uj (12) 

A 2 = Z- 1 dujp(uj) (oj — E) 2 . (13) 

Typically the largest eigenvalues reveal the ground state. 
There are few exceptions. The spectra at r = la and r = 
8a, shown in Fig. [51 are examples. Selecting the support 
5 n of a peak can be slightly ambiguous, particularly if 
the peak sits on background noise, or if there is overlap 
with another one. We have therefore, as a matter of 
course, performed least-square fits to the average spectral 
functions p(u>) with one or two Gaussians, as required in 
each case. These fits then give values for Z, E, and A 
as defined in (lll |) -(fT3 |) directly from the parameters of 
the Gaussians. In this manner we have obtained 16 sets, 
one for each r, of five-quark hadron ground state energies. 
Each set contains four energies corresponding to hopping 
parameter values k = 0.140, 0.136, 0.132, 0.128. 

The spectral density functions p(ui) typically exhibit 
some additional structure at the high end of the lu range, 
a broad peak in most cases. These secondary peaks ap- 
pear because diagonalizing the correlator matrix sepa- 
rately on all time slices only ensures that its eigenval- 
ues describe a sin gle state from the physical spectrum at 
large values for t [14J. At small time slices, close to the 
source, the eigenvectors of C(t,to) may still describe a 
mix of spectral states. In practice the large-t rule means 
that, in each eigen channel, one should only take the low- 
est peak into account. In physical terms the secondary 
peaks are separated from those by «3GeV, and larger, 
and thus should probably be considered lattice artifacts. 

In order to relate our results to the physical hadron 
spectrum we have also computed 3x3 correlation matri- 
ces with standard local Q pseudoscalar meson (it), vec- 
tor meson (p), and nucleon (N) operators, using the same 
smearing and fuzzing prescription. The analysis was per- 
formed with the MEM in exactly the same way as de- 
scribed above. This now allows extrapolations of hadron 
masses to the physical pion mass region (am^ « 0). We 
use the extrapolation model introduced in [111 ] 

y = ci + c 2 x + c 3 ln(l + x) , x = (am^) 2 , (14) 

and y = am. For a motivation of the same we refer 
the reader to said reference. The extrapolation prescrip- 
tion is, of course, a source of systematic error on our 
results. In Fig. [7] we show the vector meson and nucleon 
masses as a function of the squared pion mass and the 
extrapolations according to (fl4|) . The extrapolated val- 
ues am p = 0.293 and am^ = 0.531, as x — > 0, are used 
to set the reduced mass of the pN system to its experi- 
mental value rripN = 424.7MeV. This results in a lattice 




0.0 1 1 1 1 1 1 1 1 1 1 1 1 1 

0.0 0.2 0.4 0.6 0.8 1.0 1.2 
(am ) z 

FIG. 7: Extrapolations of the vector meson and nucleon 
masses y = am based on the model (I14[) and four hopping 
parameter values, top panel. The bottom panel shows two 
typical examples of five-quark hadron ground state energies 
y — aE, at distances r = la and r = 4a, and their respective 
extrapolation using the same procedure. 



constant of a" 1 = 2251MeV or a = 0.088fm. The cor- 
responding vector meson and nucleon masses come out 
as 659McV and 1196MeV, respectively, which deviate by 
15% and 27% from their experimental values. We take 
these numbers as indicators of the systematic errors to 
be expected due to the extrapolation. 

Also in Fig. [7] two representative examples of extrapo- 
lations of five-quark hadron ground state masses y = aE 
are shown, the relative distances are r = la and r = 4a. 
The error bars represent spectral peak widths A accord- 
ing to (fl~3)) . We utilize those to compute uncertainties 
on the extrapolated energies: At each of the four data 
points stochastically independent Gaussian random num- 
bers with average and standard deviation given by the 
energies and peak widths of the data, respectively, were 
generated and again fitted with (fT4)) . Repeating this a 
large number of times then gives rise to a standard de- 
viation which we interpret as a MEM peak extrapolated 
width A. The extrapolated energy spectrum obtained in 
this way is displayed in Fig. [5] as a function of the relative 
distance r. 
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FIG. 8: Extrapolated ground-state energy E of the heavy- 
light K—A like system versus the relative distance r with MEM 
based uncertainties (plot symbols). The solid line represents 
a fit with the model (|16|l . The dashed line shows the same 
fit, but the mirror term aV(L — x) has been dropped. 



TABLE I: Comparison of the statistical errors S of the 
ground-state peak positions E, see (|12p , and the peak widths 
A, see (|13[l. Results are shown for the lightest quark mass, 
k = 0.140, and for all relative distances r. 



la 


aE 


aA 


o S 


1 


1.193 


0.050 


0.038 


2 


1.434 


0.106 


0.031 


3 


1.514 


0.083 


0.037 


4 


1.526 


0.062 


0.059 


5 


1.512 


0.067 


0.076 


g 


1.378 


0.053 


0.032 


7 


1.206 


0.052 


0.024 


8 


1.120 


0.035 


0.029 


9 


1.190 


0.055 


0.030 


10 


1.223 


0.065 


0.020 


11 


1.348 


0.066 


0.034 


12 


1.390 


0.072 


0.030 


13 


1.465 


0.047 


0.024 


14 


1.338 


0.080 


0.030 


15 


1.315 


0.058 


0.040 


16 


1.297 


0.062 


0.009 



The energies shown in Fig. [5] are also subject to statis- 
tical errors which originate from the gauge configuration 
ensemble. To identify those a bootstrap error analysis 
[l7| was performed. Using a resampling factor of eight 
on the 343 gauge configurations the MEM s pec tral anal- 
ysis was repeated for each bootstrap sample [29[. In each 
case the same fixed random number sequence was em- 
ployed in order to eliminate fluctuations in the annealing 
process. The statistical errors 5 thus obtained are com- 
pared to the peak widths A in Tab. Q] Results are for 
k = 0.140 (lightest quark mass). The other n values give 
very similar numbers because the underlying set of gauge 
configurations is the same. Clearly the statistical errors 5 
are consistently smaller than the peak widths A. On av- 
erage over all relative distances r the ratio A/5 turns out 
to be slightly larger than two. Thus, judiciously keeping 
the larger errors, we have used the spectral peak widths 
as the principle input for computing uncertainties. 



Thus lattice sites with relative distances x and Lj, — x 
are equivalent. However, as the two-hadron operators for 
the K- A like system are essentially loop operators, see 
Fig. [21 those do not share the periodic behavior of the 
fundamental lattice fields. This effect should be most 
notable for distances around « L%/2 — 16. Second, in- 
teractions with hadron images in adjacent copies of the 
lattice are present. Their dominant effect is a constant 
interaction energy, as the relative distances of most of 
the mirror images remain constant with varying x. Thus 
motivated, empirically one finds that the lattice data are 
well represented by the model 

oVl{x) = o>q + aV(x) + aV(L — x) , x — r/a , (16) 

where «o and L are additional parameters. The resulting 
fit is shown in Fig. [8] as a solid line. The dashed line rep- 
resents q;o + aV(x) only. Table HT1 contains the fit param- 
eters. Note that L deviates somewhat from 32, as should 



ADIABATIC POTENTIAL 

The data points in Fig. [8] connect smoothly along a 
continuous curve. In order to find an appropriate fit 
model we note that the ground-state K-A like system 
should be in a relative S-wave. Our choice thus are S- 



wave 3d harmonic oscillator functions [18j. The first few 
terms are given by 

aV(x) — exp(— a\x 2 )(a2 + a^x 2 + a4X 4 ) , x — r/a , 

(15) 

with parameters oi\...a^. This will not be sufficient 
though for two reasons. First, periodic boundary condi- 
tions have been imposed on the lattice in the 3-direction. 



TABLE II: Results for the fit parameters using the model 



ao 

Ql 
a 2 

Oi3 

L 



1.213(93) 
0.0409(37) 

-0.89(10) 
0.109(19) 

-0.00328(52) 

26.54(57) 



be expected for the reason stated above. The x 2 P er de- 
gree of freedom is 0.63, the uncertainties given stem from 
the covariance matrix of the Levenberg-Marquardt algo- 
rithm. We list those for completeness only, they will not 
be used in the subsequent analysis. Rather, because we 
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are mostly interested in uncertainties of possible bound 
state energies, as the case may be, we proceed as follows: 
Each data point in Fig. [5] is randomized by stochastically 
independent Gaussian random numbers with average and 
standard deviation given by the energies and peak widths 
of the data, respectively. Then the fit with the model (flU)) 
is repeated, the corresponding set of parameters gives rise 
to a new potential V(x). The procedure is performed 32 
times. We show the results in Fig. [9] along with the po- 
tential obtained directly from the original data points. 
Specifically, Fig. [5] displays adiabatic potentials 

V a (r) = V(r/a) (17) 

on physical scales given by the lattice constant a, as de- 
termined above. 

0.6 




-2.4 ^ 1 1 1 1 

0.0 0.5 1.0 1.5 2.0 

r[fm] 



0.6 




0.0 0.5 1.0 1.5 2.0 

r[fm] 



FIG. 9: Adiabatic potential V a (r) resulting from the lattice 
simulation in physical units, panel (a). Panel (b) shows po- 
tentials from 32 sets of randomizations of the ground state 
energies E as seen in Fig. [8] 



BOUND STATES AND WAVE FUNCTIONS 

Figure [9] represents the end result of the lattice simula- 
tion. Although we expect that systematic errors coming 
from extrapolations may be as large as 30%, and the un- 
certainties on the potential itself according to Fig. b) 
are sizable, we believe that important physics can still 
be learned on a qualitative level. First, we observe that 
V a {r) of Fig. OJa) has no repulsive regions. This is a 
reasonable, and confidence inspiring, outcome because 
repulsion in a two-hadron system is typically Pauli re- 
pulsion, but of course the latter is absent in a K-A like 
system. 

Also, as mentioned earlier just before Eqn. ©, the 
relative distance r = was excluded from the simulation 
because "color twisted" operators would then have to be 
included for all distances r, to be consistent. Without 
further work it remains an open question how such a 
family of operators might influence the form of V a (r). 
However, it seems reasonable to assume that those will 
not play a major role for large r because there they would 
create hadrons inconsistent with confinement, whereas at 
r = mixing with those states would lower the ground 
state energy according to general principles of quantum 
mechanics. Thus, we believe that the current simulation 
correctly captures the most important physics of the four- 
quark system. 

The salient features of V a (r) are the two distinct re- 
gions of attraction around r ~ 0.5-1.0fm and r < 0.2fm. 
An obvious question is whether the attraction is suffi- 
cient to generate bound states of the five-quark molecule. 
To find out we have, first, computed scattering phase 
shifts 5i(p) employing the adiabatic potential V a (r) in a 
Schrodinger equation. For the reduced mass m we chose, 
in turn, the three values of the physical K-A, D-A c , and 
B-Ab systems Q. These evaluate to am = 0.152,0.415, 
and 1.172, respectively. Continuum boundary conditions 
were implemented by solving a Volterra integral equation 
within the Jost function formalism [l9| . The resulting 
scattering phase shifts are shown in Fig. |TD] Accord- 
ing to Levinson's theorem, see [l9T ] for example, we have 
5e(0) — Si(oo) — wr where n is the number of bound 
states. In the S-wave n is one for the K-A system, two 
for the D-A c system, and four for the B-Ab system. We 
need to caution though that relativistic effect are large 
in K-A, with increasing quark mass (s — > c — > b) this 
becomes less of a concern. The scattering phase shifts 
arc quite featureless, but they give us an exact count of 
the number of bound states present in the simulation. 

In order to compute bound state energies and wave 
functions we solve the radial Schrodinger equation with 
V a (r). Our conventions can be gleaned from 

-^ + [2mV a (r) + ^^-]u a = 2mEu a , (18) 
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FIG. 10: Scattering phase shifts (S-wave) for the K-A, D- 
A c , and B-Ab systems (bottom to top, respectively) versus 
the relative momentum p. 



and 

/>oc 

/ dr\u a (r)\ 2 = 1. (19) 
Jo 

The region between r = and r — 30a was subdivided 
into 60 discretization intervals. A suitable cubic B-splinc 
representation [2fj| for u a (r) was employed to enforce 
bound-state boundary conditions and provide smooth in- 
terpolation between the knots. This then leads to a gen- 
eralized matrix eigenvalue problem. The scale neutral 
bound state energies are compiled in Tab. IIIII for the 
three reduced masses of K-A, D-A c , B-Ab, and three 
partial waves. The errors come from repeating the cal- 
culation with 32 randomizations of the potential, as seen 
in Fig.H». 

TABLE III: Scale neutral bound state energies aE n of the 
adiabatic potential V a (r) for three values am of the reduced 
mass m and partial waves I. The reduced masses correspond 
to K-A, D-A c , and B-Ab- 



i 


am 


aE 1 


aE 2 


aE-i aEi 





0.152 
0.415 
1.172 


-0.22(12) 
-0.32(13) 
-0.39(13) 


-0.04(13) 
-0.26(16) 


-0.15(12) -0.001(88) 


1 


0.152 


-0.11(11) 








0.415 
1.172 


-0.28(13) 
-0.38(13) 


-0.14(12) 




2 


0.152 










0.415 
1.172 


-0.20(13) 
-0.35(13) 


-0.11(12) 





The S-wave bound state wave functions u a (r) and their 
squares |M Q (r)| 2 are displayed in Figs. 1 1 1 1 121 and fl3l The 
wave function of the K-A system, Fig. [TTJ peaks at 
around r ~ 0.6-0.9fm. This is a typical nuclear physics 
distance scale, and clearly points to the structure of a 



6 





aE = 0.22 









0.0 0.5 1.0 1.5 2.0 

r[fm ] 

FIG. 11: S-wave bound state wave function \ = u (r) (dashed 
line) and its square \ ~ \ u a(r)\ 2 (solid line) for the K-A 
system. 

hadronic molecule. With increasing (heavy) quark mass 
the number of bound states also increases. The ground 
state of the D-A c system, Fig. [TJJ exhibits the same fea- 
ture. However, the excited state wave function has siz- 
able support at small distances r < 0.2fm. The wave 
function starts to 'feel' the short-range attraction of the 
potential V a (r), see Fig. [9ja). This property is fully de- 
veloped in B-Ab, see Fig. [131 The ground state wave 
function 'lives' in the large-distance trough of V a (r) and 
thus points to a hadronic molecule. On the other hand, 
the wave function of first excited state 'lives' in the small- 
distance regime of V a (r), which indicates a very tight 
spacial structure. There is nothing as distinct about the 
third state, and the forth one strongly resembles a con- 
tinuum state. 

The wave functions of the ground and first excited B- 
Ab states with aEi = —0.39 and aE 2 = —0.26, respec- 
tively, shown in Fig. [13] deserve further comment. View- 
ing the system as consisting of five quarks, the ground 
state is clearly a pure two-hadron molecule, and the rms- 
radius evaluates to r rms sa 0.74fm. The excited state, on 
the other hand, is spatially compact in a dramatic way: 
Most of the probability density |u a (r)| 2 is within a region 
of less than 0.3fm radius and peaks just above O.lfm. In- 
evitably, this five-quark state should be interpreted as a 
pentaquark. 

We should remark here that our original five-quark op- 
erators ([I]) are indeed capable of describing such physics, 
because the set comprises operators where the heavy 
quark and anti-quark have very small distances, the 
smallest one being the lattice constant a — 0.088fm. 
Also, the B-Af, results are probably our most reliable 
ones, because relativistic effects are not very important. 

On the other hand, systematic errors, having to do 
with the adiabatic approximation, extrapolations to zero 
pion mass, and dismissal of relativistic effects, clearly ren- 
der our interpretation a qualitative one. It is interesting 
to ask though if the reason for our inability to establish 
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FIG. 12: S-wave bound state wave functions \ = u a (r) 
(dashed lines) and their squares \ = l w a( r )| 2 (solid lines) 
for the D-A. c system. 



a clear experimental signal for a pentaquark [2JJ, |22| , and 
also to irrefutably establish a lattice QCD pentaquark 
23L 24], HH , could be rooted in the assumption that the 
pentaquark shall be a ground state of bound light quarks. 
There is no a priori reason for this. The well-publicized 
quark model decouplet particle Q may well have a 
molecule-like structure. 

Thus we speculate that the pentaquark, if we under- 
stand it as a tight five-quark bound state, may not exist 
for all light flavored (u,d,s) quarks, but there is still hope 
in the heavy-flavor sector. However, a pentaquark may 
then reveal itself only as an excited state, the ground 
state being a hadronic molecule. 



Finally, as a matter of course, we show in Figs. fT4l and 
[15] the P- and .D-state wave functions of the calculation. 
And in Tab. IIVI we have compiled physical mass spectra 
of the considered five-quark systems by subtracting the 
computed binding energies from the experimental single 
meson and baryon masses Q. We give those numbers 
for completeness only, without claiming quantitative rel- 
evance. 



X 2 - 



-2 



0.0 



X 2 



0.0 



X 2 



-2 



0.0 



1.0 
r[fm ] 



: 



2.0 



aE = 0.15 




0.5 



1.0 
r[fm 



2.0 




iE =-0.26 



0.5 1.0 1.5 

r[f m] 



2.0 





aE 1= -0.39 






I I I 1 , , , 



0.5 



1.0 
r[fm] 



1.5 



2.0 



FIG. 13: S-wave bound state wave functions \ = u a(r) 
(dashed lines) and their squares x = I'MaMj 2 (solid lines) 
for the B-Ab system. 
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FIG. 14: P-wave bound state wave functions \ = u a (r) 
(dashed lines) and their squares \ ~ \ u a{ r )\ 2 (solid lines). 
The insets for am and aE n match the entries of Tab. IIIII 



FIG. 15: D-wave bound state wave functions \ = u a {r) 
(dashed lines) and their squares x = l u a( r )| 2 (solid lines). 
The insets for am and aE n match the entries of Tab. IIIII 



SUMMARY AND CONCLUSION 

We have performed a lattice simulation of a five-quark 
K-A like hadronic system on an anisotropic and asym- 
metric lattice. The heavy quark anti-quark pair was 
treated as static. Thus it was possible to compute the 
total energy of the system as a function of the relative 
distance r between the hadrons. The maximum entropy 
method was employed toward this end. An extrapolation 
into the physical pion mass region was also performed. 
The objective was to extract an adiabatic potential V a {r) 
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TABLE IV: Physical five-quark hadron masses m n , in GeV, 
for three values of the reduced mass m corresponding to the 
systems K-A, D-A c , and B-Ab. The layout matches that of 

Tab. urn 



t 


m 


mi 


m 2 


m 3 


"'4 





0.342 
0.934 
2.639 


1.11(28) 
3.44(29) 
10.02(30) 


4.07(29) 
10.33(35) 


10.58(27) 


10.90(20) 


1 


0.342 
0.934 
2.639 


1.36(26) 
3.54(29) 
10.05(30) 


10.60(27) 






2 


0.342 
0.934 
2.639 


3.70(29) 
11.11(30) 


10.66(26) 







for the relative motion. The potential turns out to have 
two distinct attractive troughs, one at intermediate dis- 
tances, r « 0.7fm, and one at short range, r < 0.2fm. 

To study the dynamics of the two-hadron system we 
have used V a (r) in a Schrodinger equation for three val- 
ues of the physical reduced mass, corresponding to K- 
A, D-A c , and B—Ab- With increasing heavy-quark mass 
m s < m c < m& the number of S-wave bound states in- 
creases from one to four. We have examined the corre- 
sponding wave functions to study the nature of the five- 
quark systems. 

Systematic errors predominantly originating from us- 
ing the adiabatic approximation, extrapolation to zero 
pion mass, the non-relativistic framework, and possi- 
bly the quenched approximation, render the results of 
our study qualitative only. Nevertheless, we believe 
that the results for the B-Ab ground and excited state 
wave functions, respectively, bring to light a particu- 
larly interesting scenario with regard to five-quark system 
physics: The ground state is best described as a hadronic 
molecule, with a relative hadron-hadron distance match- 
ing a nuclear physics scale, whereas the excited state ex- 
hibits a wave function with support on very short dis- 
tances of r < 0.2fm, or so. This leads us to interpret 
the excited state as a pentaquark. The results for the b- 
quark system are less prone to relativistic corrections and 
therefore are our most reliable ones. The light-quark K- 
A exhibits one ground state which clearly is a hadronic 
molecule. 

In the light of our results it might be worthwhile to 
initiate future lattice QCD studies with this scenario in 
mind, but with a more realistic set of operators, lattices, 
actions, and lighter quark masses, etc. All things con- 
sidered, there is ground for the hypothesis that an all 
light-quark pentaquark may not exist, but that the exis- 
tence of a genuine pentaquark in the heavy-quark sector 
cannot be ruled out. However, if present, it may well be 
an excited state. 

This material is based upon work supported by the 
National Science Foundation under Grant No. 0300065. 
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